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ABSTRACT 

This report describes the work carried out under Contract NAS8-29945, 
"Conversion of HILTOP and ASTOP Computer Programs, " The study entailed 
the detailed documentation of the programs, the delivery of IBM 360 versions 
of the programs to the NASA Marshall Space Flight Center, and the investigation 
of selected topics relating to the possible extension of the HILTOP program. 

The documentation of the computer programs is being published concurrently 
with this summary report, and the programs have been delivered. This report 
presents the analyses of the possible HILTOP extension and also gives the results 
of an extra-ecliptic mission study performed with HILTOP. 
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I. INTRODUCTION AND SUMMARY 


The primary purpose of this study was to develop detailed documentation 
for two existing low thrust trajectory optimization computer programs: (1) HILTOP, 
a variational calculus program for optimizing heliocentric trajectories in a central 
force gravitational field; and (2) ASTOP, a parameter optimization program for 
n-body interplanetary trajectories. 

Prior to documenting HILTOP, a concerted effort was undertaken to reduce 
the core requirements of the program to facilitate its use at installations where core 
availability is at a premium. This reduction in program size was achieved through 
various techniques, such as reduction of array dimensions and elimination of 
selected subroutines. None of the modifications made restrict or degrade the pro- 
gram's capabilities, although a small degradation in efficiency will result for some 
mission types. For example, separate routines to handle two-dimensional missions 
were removed. The core requirements were reduced about 25 percent. The docu- 
mentation ^ is being published concurrently with this summary report. 

Much more extensive changes were made to ASTOP before and during docu- 
mentation. Unlike HILTOP, ASTOP is a research program which had been used 
very little and contained much unused code that remained from the ITEM program 
from which ASTOP was derived. The first step undertaken was to pass the source 
code through various utility programs which reordered the non-executable state- 
ments (such as type statements, dimensions, commons, data and equivalence 
statements), sequenced the statement numbers and, in general, cleaned-up the code. 
This was followed by passing the new source code through additional utility programs 
which generated various cross reference tables of subroutines, labelled commons 
and common variables. A careful scrutiny of the compilation listings in conjunction 
with these cross reference tables resulted in the elimination of certain subroutines, 
commons, blocks of code and selected variables from the program because they 
were not used. The program was then debugged and forced to execute a test case 
to assure that no serious oversights had occurred. The entire process was then 
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repeated, again resulting in a significant reduction in non-productive code, and the 
documentation of subroutines was then begun. During the documentation process, 
virtually every line of code was studied and several additional changes were effected. 
In all, the region size was reduced from about 480 K to 336 K. The documentation 
for ASTOP is given in Reference [ 2 ] . 

The study also called for the investigation of several potential extensions to 
the HILTOP program. These were: 

(1) extension of fixed thrust cone angle logic to include multiple cone angles; 

(2) generalization of fixed thrust cone angle logic to permit solar array 
orientation to be non-normal to the sun line; 

(3) improve the convergence characteristics of the iterator; 

(4) include solar array radiation degradation effects in the performance and 
optimization model; 

(5) generalize the propulsion system model to include variable specific impulse 
and efficiency as a function of throttling ratio; and 

(6) alleviate convergence problems manifested by infinitesimal thrust or coast 
phases. 

None of the above extensions was implemented in HILTOP as part of this study be- 
cause the magnitude of the task exceeded the funds available. Each potential ex- 
tension was studied, although at varying levels of detail. 

The items (1) and (2) above were formulated within the same general problem, 
and the detailed analysis is presented in Section II. With regard to item (3), a new 
iterator was developed; however, the approach was designed to speed the computa- 
tions on each iteration rather than improve the convergence rate. A detailed dis- 
cussion of this iterator is given in Section III. A radiation degradation model was 
actually implemented in HILTOP as part of a separate study. The intent of item 
(4) was to improve this model using inputs to be provided through MSFC. Since no 
revised model was developed, no further analysis was undertaken. A detailed 
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description of the current model is given in [l] . A model of thruster throttling 
effects was incorporated in an early version of HILTOP but, because its use was 
rather expensive, this model was not included in the current version. An analysis 
of the principal requirements of such a feature is presented in Section IV. The 
convergence problems alluded to in item (6) have been given considerable thought 
and, as yet, no satisfactory solution has been found. The present understanding 
of the problems is discussed in Section V. 

During the course of the study, an analysis was performed of extra-ecliptic 
missions using the SEP stage atop a Titan III E/Centaur launch vehicle. Several 
trajectory classes involving mission durations from I 2 to 3 years were investi- 
gated for heliocentric inclinations ranging from 45 to 60 degrees. The results 
of this analysis are presented in Section VI. 
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II. GENERALIZED ANALYSIS OF FIXED THRUST GONE ANGLES 


Optimal trajectories with unconstrained thrust direction will frequently re- 
sult in a thrust angle relative to the sun line that fluctuates over a wide range during 
the course of the trajectory. With SEP systems, for which the arrays are usually 
assumed to continuously face the sun, this requires a continual movement of the 
thrusters relative to the arrays, a requirement that is highly undesirable. For 
this reason the concept of operating the system with a fixed spacecraft array con- 
figuration is of much interest. The capability of simulating this constraint has 
been available in HILTOP for some time. However, the performance penalty in- 
curred in some missions with a single fixed cone angle is excessive, so the ability 
to define the performance sensitivity to a number of fixed angles is desired. 

The following development represents an extension to the formulation 
of the HILTOP program as presented in Reference [l] . Familiarity with 
these two documents is essential to the clear understanding of the analysis which 
follows. 

Consider the case of a solar electric spacecraft with solar array orienta- 
tion defined by the unit vector n and thrust in the direction of the unit vector 

e , and suppose that e is constrained to lie nominally at one of a number of 
t t ^ 

specified cone angles 0., i = 1, 2, , k, from ii. Also, to provide for the 

option of thrust vectoring, admit the possibility that e t may lie anywhere within 
a cone of specified half angle about the nominal directions defined by ?>.. (See 
the sketches at the end of this section). This constraint may be expressed mathe- 
matically by the inequality 

^ = (cos 1 (e t • n) - 0 . ) 2 - 7?. 2 ^ 0 . (2-1) 

In addition, it may be desirable in certain cases to orient the solar arrays to 
continuously maintain maximum power output. This may be accomplished by 
imposing the constraint 

^2 = n • e r - 1 = ° for r * r c , 
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or 

_ 2 2 

$ = n • e -r/r =0 for r < r , 

2 r c c 

where e^ = R/r and r^ Is a specified solar distance below which the arrays are 
tilted so as to maintain constant power output of the arrays for r < r^. 

To the state equations given in [l ] , we add for this problem the k equations 

0. = O, 1 = 1, 2, , k. 

These are included to yield associated adjoint variables which will appear in 
transversality conditions if it is desired to optimize the k cone angles. 

The variational Hamiltonian for this problem is written 


(A ' V C v +x r]" A ' ® 

r 

r - 2 21 — - 

+ A (cos (e *n)-0.) - 77 . + A (n • e - p) , 

x. L t 1 1 J y r 


( 2 - 2 ) 


2 2 

where p = r/r if r<r and p = l otherwise. Of course, A and/or 
c c x. 

1 

are zero if the associated constraints are not imposed. The optimal control 
problem now is to choose e , n , and h at each point along the trajectory so as 
to maximize h^ subject to the specified constraints. Since the last two terms in 
(2-2) never contribute to the magnitude of h , it is seen by inspection that h^ is 
maximized with respect to e and n by choosing e^ as close to A as possible 
and choosing n so as to make y as large as possible. Of course, any constraints 
between e and n preclude choosing e and n independently; therefore, it is, 

t l 

in general, necessary to compromise in maximizing y and (e^ • A) individually 
in favor of maximizing the function y ( A • e^ - — A^). This must be done by 
considering individual cases that may arise. 

First, consider the case for which the solar arrays orientation is con- 
strained so as to produce maximum power output. Under this constraint one 
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can consider maximizing only after the constraint is satisfied, and maximizing 

h^ is equivalent to maximizing (A* e^) subject to the constraint. Let a denote 

the angle between R and A and let j denote the index of the currently optimum 

cone angle (the determination of which cone angle is currently optimum will be 

considered subsequently). Then, for r > r , the constraint of maximum-power 

c 

output requires that n = e^, and the choice of which maximizes (A * e^), and 
therefore h , subject to the constraint is 



e cos ( 0 . + 77.) + (m x e ) sin (0 . + 77 ) if a ^ 0 . + 77 
r J T r J V J J 

e N if 0 . - 77. < a < 0. + rj. 

X J J J J 


e cos (0 . - 77 .) + (m x e ) sin (0 . - 77.) if a ^ 0 . - 7?. 
r v ) j r J ) ) ) 


where e. = A/X and m = (R x A)/|r x A J. For r < r , n is constrained to 
A c 

lie on a cone of half-angle 


n -1 , 2, 2 
0 = cos (r /r ) » 
c 


about e , and the optimal choice for e is 
r t 


e t = < 


e cos (0 . + 77. + 0) + (m xe ) sin (0 . + 77. + 0) if a s 0 . + 77 . + 0 
r J J r j 'j J j 

e s if 0 - 77. - 0 < tt < 0. + 7), + 0 (2-3) 

X 3 J J J 

e cos (0 . - 7?. - 0) + (mx e ) sin (0 . - 77 . - 0) if 0 . - 77 . - 0 
r J J r j ] ) j 


while n (which is not always unique) may be defined 


n = 


e cos 0+ (m xe ) sin 0 if a ^ 0. + 0 

r ' r 1 


e cos 6 + (mxe ) sin 0 cos € + m sin 0 sin € if 0. - 9<a<0. + 0 
r yr j J 


e cos 0- (mxe )sin0 if a ^0.-0 

r r j 


(2-4) 
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where e = cos 1 [(cos 0^ -cos 0 cosa)/sin 0 since j . Note that equations (2-3) and 
(2-4) also hold for the case r > r^ if one sets 0=0. 

For the case in which n is not constrained to continuously yield maximum 
power output of the arrays, the optimal control problem becomes one of maximiz- 
ing the function y (e^ • e - b) subject to the cone angle inequality constraint 

(2-1), where b = A, v/oX. As in the preceding case, when 0. - tj. - 6<a< 
v _ ) ) 

0. + V. + 0 the quantities y and e. • e maybe maximized independently while 
J J At 

satisfying the cone angle constraint by rotating n out of the plane of R and A . 
When a is not within this interval, both n and e t must lie in the plane of R and 
A and the maximization of y (e^* e t - b) may be taken with respect to a single 
parameter, say the angle 6 between e^ and e^. (See Figure 2 at the end of this 
section). This is accomplished by solving the equation 

c) y 

(cos 6 - b) “ - y sin 6 = 0 , 
for 6 subject to the condition 

3 3 y 

(cos6-b) — - -2 sin 6 -r-? - ycos 6^ 0 ) 

36 2 

to assure the function is maximized. The solution of the equation (2-5) for 6 
will, for most forms of y, require an iterative technique. For our purposes, 
y written in terms of 6 is of the form 


(2-5) 


cos (a— 0 -r} -6) . ((i + 4)/4) 

'■ ) a ,( r 1 " J 


i=0 


so that 


by x V / 1+4 w c °s (a-0 -?7 -6) ((i + 4)/4) 

1=0 r 


A suggested approach to the solution of equation (2-5) is to employ a Newton’s 
iteration with sin 6 as the independent variable, using as a first guess 
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sin 6 = 7-— tan (a-0. - ri.) . 

2-b j V 


Once the optimum value of sin 6 is obtained, form cos 6 = J l -sin^6 and 


write 


or 


e^ = e^ cos 6 - (mxe^) sin 6 , 


n = e, cos (0 . + 77.) - (mxej sin (0 . + 1? ) if a>0 + t?. + 0 

t J J t J j j J 

n = e, cos (0 . - 77.) - (mxe.) sin (0 . - 17.) if a<0. -17. - 0 . 

t j j t' j y j j 


Of course, if 0. - ri. - 0<a<0 +77+0, then 
J J J J 


with 




e t = V 


n = e cos 0+ (mxe ) sin 0 cos e + m sin 0 sin € , 
r r 


0 if a^0 . + 0 

J 

cos ^["(cos 0 . -cos 0coscc)/sin 0sina 1 if 0 . - 0 < Qf < 0 . + 0 4 

L- J —I J J 

IT if O!^0 . - 0 

J 


To determine which of the 0 is optimum at any instant, assume that 
0 < 0 < < 0, , and suppose that, at this instant, 

\ A K 

0. + t?. + 0 < Ct < 0 . „ -t). -0. 

1 '1 i +1 'i+l 

Then j, the index of the optimum cone angle at that instant is 

f 1 if ( 0 i+i“ Vi“°° ~ ( a - 0 i-V >o 

3=1 


m lf ( Vi-Vi' a) ' ( 0 ! "V , V <0 
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The switch from one cone angle to the other occurs when the difference vanishes. 


Since h is linear in h , the choice of h is made as described in 
v <7 a 

[ 1 ]. 

The adjoint equations are obtained formally through partial differentiation 
of the variational Hamiltonian. Those that differ from the adjoint equations pre- 
sented in [3] are 


A = 2M (A-E)R- (A • e - - A ) + h X 1 

5 ' ' 3 r L a V t c V ’ p yj 


[S-s(S • 5> 5 ] 


where 


and 


X" . = 2X [cos 1 (e, - n) - 0 . 1 ; X* = 0 
0. X. L t 3 J 0. 

J J 1 

1 X / >^ \ / e * n \ */ 4 
r 1=0 r 

h ={° ifr>r c . 

p 1 1 if r < r^ 


for i ^ j , 


The Lagrange multipliers X and X are determined by setting the 

X i ^ 

variations in h resulting from independent variations in e and n, respectively, 

V t 


to zero. That is, 

Ih ^A-2A 
L a v 3 


( cos 1 (e * n) - 0 . ) 

■ - v ■, ] "]• a v°- 

i Jl - {e t • n) 


{ fli_ (A’ e,--X ) + X ]e -2X 
i L a v t c v yj r x. 


(cos 1 (e t • n) - 0 . ) 


«A-<V” )Z 


J- 6 " = 0 - 


(2-6) 


(2-7) 


Now, because the variation of a unit vector must be normal to the unit vector, 
it is clear that 6e may be divided into two components - one along (n x e ) 

t t 
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and the other along (n x e ) x e . Since variations in these two directions are 
independent, the equation containing 6e^ must be satisfied by the variation 
along each component independently. Substituting into (2-6) the variation along 
nxe^ and using the identity n • (n x e^) = 0 leads to the result 

A • (ii x e t ) = 0 , 

which indicates that A, n and e^ are coplanar vectors. Then, substituting 

into (2-6) the variation along the second component yields the desired definition 

of X . Employing the identities 
x i 

A • [(n x e t ) xej=- (n x e fc ) - (A x e t ) f 
n * [(n x e t ) x e t ] = - (n x e t ) • (nxe^-fl-^ 1 ii) 2 ) , 
m = (n x ej / | n x e^ | » 


yields for 




gr 

v 


m* (Ax e ) 

1 / 


2 [cos 1 (e • n) - 0 . ] 
t j 



for i ? j . 


Note that the identity involving m is valid only outside the interval 


V 


77. - 9 < a < 
) 


. + n +6 
j j 


but, since A x e^ becomes the null vector when a is within the interval, the 
above expression for X is valid for all a. 

j 

Before defining X recall that X is non-zero only if the array 
orientation is constrained to yield maximum power output. Also note that X 
appears only in the equation for A where it is multiplied by the step function 
h . Therefore, X^ influences the problem only when the array orientation 
constraint is imposed and when r < r , and we will confine the discussion of 
X to cases where those conditions apply. Proceeding as with 6e^, consider 
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fin broken down into the two components along (n x and [(n x e t ) x n], Em- 
ploying in (2-7) first the component along (n x e ) and noting the identity 

l 

e * (n x e ) = 0, one is left with the condition 
t t 

fh (A * e, - - X ) + X I fe 5 (n x e, ) 1 = 0 . 

L a v t c v yJL r ' t J 

Now, when a is outside the interval 

0.-17. - 0 < a< 0. + n + 0 , 

11 11 


e , n and are coplanar such that e f * (n x e^) = 0, and no information is 
given about X^„ However, when & is within the interval, n is rotated out of 
the plane of R and e^, and X is then defined by the relation 


h (A * e - — X ) + X = 0 . 
a v ' t c v' y 


It remains to define X when oi is outside the interval, and this is done by 

y _ _ _ 

considering the component of fin along (nxe^xn. Employing the identities 

e fc * [(n x e^) x n] = (n x e^) * (n x e^) = 1 - (e t • n) 2 , 

e • C(n x e ) x n] = (n x e ) • (n x e ) , 
r t t r 

m = (n x e t )/|n x e t | , 


and substituting in (2-7) for 


X leads to the relation 
x. 
l 


h 

CT 


gy 1 

V 


(A * e - - X ) + X 
' t c v f ' 



v 


m • (A x e ) 

U 

m • (n x e ) 
\ r 


which completes the possible cases for which it is necessary to define 


X . 
y 


As a final point, it should be noted that the transversality conditions to 
be satisfied if the k fixed cone angles 0. are to be optimized are, simply. 
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\ ( V- X 0 (t o>" 0; , = 1 ' 2 - — ■ k 

i i 

where, without loss of generality, X^ (t^) maybe set to zero. 

In summary, it is helpful to note the principal modifications to HILTOP 
that are required to Implement the generalized cone angle logic. These are: 

(1) At the start of each thrust phase, determine the index j of the optimum 
cone angle as prescribed in the text above. 

(2) At each point along the trajectory, monitor the function * 


(0. + 0.^) - (V. ~ V.^) ~ 2a. 

When the function becomes positive, it is necessary to iterate to the time that 
the function is zero. At that time, decrement j by one, re-evaluate all deriva- 
tives for the new cone angle and continue integrating. 


(3) At each point along the trajectory, monitor the function* 



+ 0 .) - (77 


j + l 


- 77.) - 2a . 


When this function becomes negative, it is necessary to iterate to the time that the 
function goes to zero. At that time, increment j by one, re-evaluate all deriva- 
tives for the new cone angle and continue integrating. 


(4) Add the new terms to the derivatives A . 

(5) Add an additional equation to the integrator; i.e., that required for X^ . 

J 

Note that only one additional equation need be integrated at any instant in time 

since X* =0 for i^j. 

0. 

1 

(6) Add the additional transversality conditions, X^ (t^) = 0. 

i 


(7) Add the computation of X and X . 

X i y 

(8) For cases in which n is not constrained to lie along e^, the 
solution of Equation (2-5) is required at every derivative evaluation. 


* Du ring thrust phases only. 


iterative 
This will 
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significantly increase the CPU time required to integrate a trajectory. 
(9) Evaluate the unit vectors e^ and n. 



Fig. 1. Spacecraft geometry and nomenclature. 


ORIGINAL PA@i is 
OP POOR QUALITY 



Fig. 2. Orientation geometry and nomenclature. 
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III. GENERALIZED RANK-PRESERVING SEARCH ALGORITHM 


The following analysis was evolved by Samuel Pines and consists of a 
search algorithm designed to produce accelerated convergence in terms of signifi- 
cantly decreasing the number of operations per iteration. 

In many problems involving the solution of a system of nonlinear equations, 

A 

z, which are functions of a variable control vector, x, the required vector, x, 
is arrived at by successively reducing the residual error vector to zero in a sequence 
of linear transformations. As the solution approaches the correct solution, the 
linear transformation approximation dominates, and the required solution is ob- 
tainable in an accelerated manner. In order to ensure that all states remain reach- 
able during the iteration process, it is important that the rank of the linear trans- 
formation remains complete. 

The method outlined in this section possesses the acceleration quality, 
makes no assumptions of symmetry in the desired transformation, is rank pre- 
serving and reduces the errors in the transformation to zero in a least square 
sense. The method applies to systems of nonlinear equations in which the number 
of equations is equal to the number of control variables. 

In the analysis which follows, x is the control vector of dimension n, 

A 

z is the nonlinear n vector output as a function of x, is the desired output 

A A 

n vector, y is the residual error vector, z.(x.) - z , x is the desired control 
^ j, LIU 

vector whose output is z D (x), A is the local linear (nxn) transformation, H. 
is the approximation to the inverse of the matrix A, and oc is a positive scalar. 
Equation numbers referenced in this section pertain to this section only. 

In the derivation of the search algorithm, we assume that in the neighbor- 
hood of the desired solution, x, the problem has a linear representation in 
matrix vector form which is given by 


z D (x) = z. (x.) + A(z i (x.), x.){x - x.}. 



The matrix, A, is assumed to be unavailable, or too cumbersome to compute. 
However, we assume that we have a numerical process by which we can readily 
compute the vector output, z.(x^), for each given control vector, x^. The final 
solution is given by 


where 



y i = Z i (X i> ' Z D (i) • 


Let be an a priori, full rank, approximation to the desired A 
matrix, and let and y be the initial estimate of the control and the error 
in its output from the desired output, respectively. For any control, x we 
have the improved control, x. , given by 

L+l 


x. = x. - a. H. y. , 

i+l i i ri’ 


(3-D 


where a. is a positive numerical scalar so chosen that the resulting residual 

in the output y. (x. ), is smaller in magnitude than the magnitude of the 

i+l l+l 

previous, y.(x.). The full linear step in x. - x. corresponds to a. = 1.0. 

I L ^ l+l i l 

To obtain the next approximation in A , we produce H._^ with the following 
properties: 


{ y t + i - y t} Tl Vi T = { X L + 1 - x i) T + 6xT . 


H i+1 T = H ( W. 

2 2 

These are the n +n equations in n unknowns, and we will obtain a solution 

2 

for H. + ^ which minimizes the sum of the squares of the n +n elements of 

fix and dH.. 

i 


Let 


Ax. = x. „ - x. 
i i+l i » 


Ay i " y i+l ' y l ' 
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The least squares solution for H ls given by 


Eh, * H - - 

1+1 l 


1 7 — ( H i A V Ax [} Ay t' 


1+ Ay Ay 


To investigate the rank of H. + ^, we note from Eq, (3-1) that 

H i + i = H i (’ - , S: - Wi - ay J Ay i T )' 

1 + Ay Ay 


Thus rank (H ) = rank (H.) provided the determinant 


det ( I 77— { Ay i - ay J Ay i T ) * 0 

1 + Ay Ay 


(3-2) 


The determinant will not vanish provided 


a ^ 


A T 

Ay i y i 


But, if for a positive a, |y | < |y. | we have 


T T 

cx. > 0 > y. , y. - y. y. > 

J i+1 J i 


so that the rank preserving inequality of Eq. (3-2) is always satisfied for an 
acceptable correction step in the control. 

The advantage of this iterator is that once a reasonable estimate of an 
initial H is available, no matrix inversions are required nor is there any 
necessity to generate, either through analytical equations or differential 
corrections, a matrix of partials of the dependent parameters with respect to the 
independent parameters, A limitation of the iterator is that the number of de- 
pendent and independent parameters must be equal. In any specific application, 
the most critical factor, in terms of overall efficiency, will be the algorithm for 
controlling the parameter a. 


17 



This search algorithm has been incorporated into a breadboard computer 
program using the identity matrix as an initial estimate of H and a simple control 
law for a, and its convergence ability has been demonstrated,, The logical sequence 
of the experimental computer program is given in the flow chart presented on the 
following page. The full potential of this new search algorithm has yet to be ex- 
plored . 
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Flow Diagram 

Rank Preserving 
Search Algorithm 





I lyj ! = I |y 2 l I 


P = (l + Ay T Ay) -1 

to= p (HAy - Ax) 

H = H - 4oAy T 
v = Hy. 
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IV. TRAJECTORY OPTIMIZATION CONSIDERING THRUSTER THROTTLING 


With SEP, the power generated by the solar arrays and input to the power 
conditioning units varies continuously with distance from the sun. The utilization 
of this varying power by the thruster subsystem will, to some extent, vary with 
the assumed design of the thruster subsystem. The model incorporated in the 
HILTOP program assumes constant efficiency and jet exhaust speed throughout 
a mission. Since hardware design considerations preclude these operating condi- 
tions, except possibly at non-optimal settings, it is desired to learn the effects on 
performance of simulating a more x-ealistic thruster subsystem control policy. 

The formulation of such a simulation capability in the HILTOP program is pre- 
sented in this section. Again, knowledge of the nomenclature and formulation of 
HILTOP as presented in Cl] is assumed. 


This analysis assumes that the thruster subsystem is comprised of a number 
of individual thrusters, each of which is characterized by a reference power p^ 
which represents a maximum allowable input power. The throttling ratio q is 
defined as the power input to the thruster subsystem divided by the total allowable 
power input to the n^ operating thrusters, i.e., 


p ref y 

Vt 


where p is the reference power available to the power conditioning units at 
ref 

1 AU from the sun and y is the ratio of power at a distance r to power at 1 AU. 
This implies that all thrusters currently in use are operating at the same throttling 
ratio. The formulation employed is based on the assumption that individual thrusters 
are turned on or off as necessary to make optimum use of the power that is available 
and to satisfy any constraint on a minimum allowable throttling ratio q , 

The variations in jet exhaust speed c and efficiency r) are assumed to be 
specified functions of q. The precise form is not important here; we will simply 
write 
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c = c(q); T 7 = 7?(q) , 
with the boundary conditions 

^=0(1); 7? 1 = n(i). 


applying at q=l. 

[ll 

Employing the nomenclature of the HILTOP users’ manual , the formula- 
tion for the throttling capability is as follows. Define at each point along the tra- 
jectory the number of operating thrusters as the smallest integer such that 


■t* 


’’ref 7 ' 


Then the throttling ratio of each operating thruster is defined 


q = 


P ref y 

Vt 


If, however, this value of q is less than the input q , it is necessary to de- 
crement n^ by one, after which 


n p < 
t*t 


P ref y ’ 


and q is reset to unity. This results in the utilization of less power than is 
available to the power conditioning units, and it is assumed that this excess power 
is radiated off in space* with no penalty in performance. Then the thrust accelera- 
tion that the propulsion system is capable of generating at any instant may be 
written 

2n t p t q7, 

a = — . 

m c v 
o 

For comparison, the thrust acceleration in the HILTOP formulation was written 

a y 
o 

V ° 

*Or the solar panels are tilted so that only the power required is produced. 
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The optimal switching of individual thrusters is determined from the maximum 
principal. The variational Hamiltonian is written 

V h cr[ a(A ’ Vo A ’ R -A-R. 

r 


Defining the switch function o as 


where 


<7 = cr* + A^/a , 
a* = A • e t - v\^/c , 


it is seen that 

fO if a < 0 
a 1 l if a > 0 * 


When h = 1, optimal switching of individual thrusters is governed by choosing 
the number that maximizes a a*. If y < 0, clearly the choice is whether to 
continue using the current number of thrusters or to turn one off. Conversely, 
if y > 0, the question is whether another thruster should be turned on. Con- 
sider first the case where y < 0 in which one looks for the condition 


n t qr ? 

c 



( V l) 




*7, 


9 


to be satisfied. When this occurs, n^ is decremented and the thrust acceleration 
is discontinuous although the Hamiltonian and all adjoint variables are continuous. 
After decrementing n^_, the condition , 


n t p t < 


p ref y> 


will exist for a finite period of time during which the power utilized is held con- 
stant, and the power differential p y-n p is radiated into space. This con- 

tGl t L 

dition is maintained until y decreases to the point for which the inequality becomes 
an equality after which the thrusters are again throttled as necessary. 
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o 

The case for y > 0 is basically the reverse of the above. The n^ thrusters 
are throttled up until q = 1 is reached. Thereafter constant power is utilized until 
the condition 


Vi 




is satisfied where 


p ref y 


% ' (n t + l)P t • 


and t] and c are the efficiency and jet exhaust speed, respective^, evaluated 
P P 

for q = q p _. 


The above conditions for optimally switching individual thrusters apply if, 

• • 

at the switch point, q for y < 0 or q for y>0 exceeds q . If this is not 

p m 

the case, then the switch is made when the condition 


q=q m if y<0 • 
or 

V q m “’' >0 ’ 

is satisfied. When this occurs, the quantity acr* is discontinuous, but the 

Hamiltonian remains continuous due to a corresponding jump in the primer 

derivative A. Letting the superscripts - and + denote limiting conditions before 

* 

and after the switch, respectively, the discontinuity in A may be written 

R • R 


This equation applies also when the last thruster is switched on or off, at which 

^ 4 " — 

time either a or a is zero. Note that the decision to switch on or off indivi- 
dual thrusters must be made after every integration step, and if the decision is 
to switch, an iteration to the switch point is required prior to continuing integration. 
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Due to the dependence of c and rj on q which, in turn, is a function of 
r, the second order differential equations for the primer vector become 



2p ref’ ?7 ' 

m ^cr 
o 



1_ dc 
c dq 



q 

_ v 

2 

c 


dc 1 

dq J 


R 


+ ^ (A" R)R - ^ A. 

r r 

Furthermore, since c is strictly a function of q, there is no point in including 

c as a state variable. Therefore, the differential equation for X may be 

c 

eliminated. 
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V. INFINITESIMAL TKRUST/COAST PHASES 


In the performance of a typical mission study, one frequently desires to 
generate a sequence of optimal trajectories over a range of values of one or 
more independent parameters. A standard approach in such a study is to use the 
results of one converged case as initial guesses to the solution of the next case in 
the sequence. When employing a variational calculus program, such as HILTOP, 
in this situation, a convergence problem occasionally arises that deserves special 
attention. If, as one progresses through the sequence, a new thrust or coast 
phase appears in or disappears from the solutions, it is not uncommon to experience 
severe difficulty, or failure, in convergence near the transition solution. Although 
this problem will arise in only a very small percentage of cases, the analyst will 
find that such problems will consume the greatest proportion of his and the com- 
puter's time over the course of the data generation task. Consequently, the 
solution of this problem would yield a significant cost reduction in electric pro- 
pulsion mission analysis. Unfortunately, no real progress has been made in 
solving this problem. The following paragraphs explain the behavior of HILTOP 
as the problem is encountered. 

Consider, as an example, a sequence of cases in which the switch function 
passes through a positive minimum in the same general region of the several cases 
in the sequence. Suppose that, as one progresses through the sequence, the value 
of the switch function at the minimum decreases to the point that, on the next case 
of the sequence, one would predict that the minimum value will be negative. That 
is, on the next case one would predict a small coast phase in the vicinity of the 
minimum in the switch function. 

In using HILTOP in a situation such as this, one should expect, at best, an 
increase in the number of iterations required for convergence and, at worst, 
failure to converge. The reason for this is that the behaviors of two neighboring 
trajectories, one with the coast phase and the other without, are different when 
subjected to the same perturbations in the independent parameters. Another way 
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of saying this is that the partial derivative matrix wili not be continuous or smooth 
across the transition. 

To understand the behavior of HILTOP when this situation arises, it is 

necessary to consider the characteristics of the iterator MINMX3. Given a vector 

AY of desired changes in the end conditions, the partial derivative matrix P, and 

two arbitrarily defined diagonal positive -definite weighting matrices W and W , 

x y 

the independent parameter correction vector Ax is obtained by solving the follow- 
ing set of simultaneous equations. 

(P T w p + Aw ) Ax- p T w Ay 

y X y 

The scalar X is termed the inhibitor and it controls the step size taken on a 
single iteration. Note that if X = 0, the formula degenerates to the standard 
Newton-Raphson algorithm. Conversely, as X approaches infinity, the correction 
AX approaches the null vector since the right hand side is finite. On each itera- 
tion MINMX3 automatically adjusts X to control the step size on the next iteration. 
For this purpose, the scalar quantity 

q = Ay t w Ay, 

y 

which represents a weighted sum square of the residuals, is formed. The iterator 
is designed to require that no trial trajectory may be accepted as the next nominal 
unless the value of q for that trial is less than the value of q on the last nominal. 
Since one is guaranteed an improvement in q by setting X to a sufficiently large 
value (i.e., by taking a sufficiently small step), the step size control algorithm 
consists simply of the logic required to successively increase X until a value of 
q is achieved that is less than that of the nominal. Once this is accomplished, 
partial derivatives are produced for the new nominal and X is reduced by a fixed 
factor in preparation for the next iteration. 

With the above algorithm in mind, consider once again the example 
hypothesized earlier. On the last case without the coast phase, a correction 
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vector AX is evaluated which normally will introduce the coast phase. Because 
the partials are incorrect for trial solutions with the coast phase, the scalar q will 
normally be larger than for the nominal. Consequently, A is increased to reduce 
the step size. Typically, one of three situations will then develop. (1) The value of 
q on solutions with coast phases is always larger than the nominal. This results in 
the choice of nominal trajectories that are always on the same side of the transition, 
with successively smaller and smaller correction vectors resulting from larger and 
larger values of the inhibitor. Eventually, the correction vector is smaller than 
the least significant digits of the independent parameter vector and the iteration is 
terminated. (2) Nominals are selected on both sides of the transition; how’ever, 
partial derivatives evaluated on both sides drive the solution toward the transition 
case where the switch function passes through a minimum with a value nearly zero. 
This condition terminates the iteration in the same manner as above. Conceptually 
in this case the problem is caused by the function q not being a smooth, convex 
function of the independent parameters. The situation in two dimensions may be 
depicted as shown in the following sketch. Note that locally the apparent correct 



choice of the independent parameter x. is in the direction of the transition solution. 
(3) A step AX is taken which yields a new nominal with a small coast phase and 
which subsequently leads to the converged solution. In the context of the hypothetical 
two-dimensional case sketched above, either of the two situations depicted in the 
sketches below would lead to the desired result. 
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To this point, no satisfactory algorithm has been found which assures that 
the third condition above is established in an arbitrary situation. The following two 
steps have been taken and are employed routinely in the HILTOP: (1) the wasting 
of some CPU time is avoided by terminating a case after an input number of tra- 
jectories have been encountered in which the magnitude of the switch function at 
a stationary point is less than a specified tolerance, and (2) consistency in the 
partial derivative matrix is maintained by adjusting the perturbation step sizes as 
necessary to ensure the same number of switch points on the nominal and all per- 
turbation trajectories. 
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VI. CURSORY TRAJECTORY ANALYSIS OF THE EXTRA -ECLIPTIC MISSION 


Due to the recent interest in the extra-ecliptic mission as a possible 
engineering test flight for SEP, a cursory investigation of launch opportunities and 
performance requirements was undertaken for this mission. The basic mission 
considered terminated in a circular heliocentric orbit of nominally 1 AU radius 
with ecliptic inclination ranging from 45 to 60 degrees. The launch year chosen 
for the study is 1979; however, the solutions obtained repeat annually and are valid 
for any launch year. 

Extra-ecliptic missions are noted for their requirement for relatively 
large geocentric declination of the hyperbolic launch asymptote. If one neglects 
the effects- of declination on launch vehicle performance, the declination will vary 
from about 50 degrees to about 70 degrees, depending on the time of year of 
launch. The lower declinations occur for launches around the vernal equinox 
(negative declinations) and the autumnal equinox (positive declinations). In HILTOP, 
the launch vehicle performance is modelled as a function of the hyperbolic excess 
speed v , the parking orbit inclination i if greater than the launch site latitude, 
and the geocentric declination 6 of the asymptote if greater than the parking orbit 
inclination. The program permits any one of these three parameters to be fixed 
or optimized independently of the other two. Generally speaking, if i is allowed 
to be optimized the resulting value will nearly equal the declination. This is be- 
cause the penalty in launch vehicle performance is much less for changing parking 
orbit inclination than for an equivalent non-coplanar injection maneuver from 
parking orbit. Certain launch vehicles are restricted to maximum parking orbit 
inclinations, however, due to range safety considerations. Such is the case for 
the Titan III E/Centaur launch vehicle which was specified for the study under dis- 
cussion here. Consequently, the geocentric parking orbit inclination was limited 
to a maximum of 36 degrees, which was used in all cases generated. The other 
two parameters, v w and 6, were optimized to yield maximum final mass. This 
approach to the inclusion of the effects of large launch asymptote declinations re- 
sults in a finite, non-zero angular offset of the initial thrust acceleration vector 
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from the launch asymptote. This offset is in the plane of the excess velocity and 
the polar axis of Earth and is in the direction of the North Pole for positive de- 
clinations and the South Pole for negative declinations. That is, the excess velocity 
always lies between the initial thrust acceleration (i. e. , the initial primer vector) 
and the equatorial plane. This offset angle is denoted Ot . 

Mission durations ranging from 1.5 to 3 years were to be investigated. 
Over such a broad range of flight times, it is essential that more than one class 
of missions be studied, A mission class for the extra-ecliptic mission may 
generally be categorized in terms of the number of burns. Typically one burn 
will occur around every nodal crossing such that, over an entire mission the 
number of burns n may be determined as a function of the number of revolutions 
r as follows: 


n = 2 r + 1 . 

Previous studies for the 1 AU circular final orbit have yielded the result that the 
preferrable class of solutions is that for which the trajectory remains in the 
vicinity of 1 AU throughout the mission. This is achieved by choosing the number 
of revolutions equal to the flight time t^ in years. For this reason mission 
classes of 4, 5, 6 and 7 burns were chosen for the flight times of 1.5, 2, 2.5 
and 3 years, respectively. 

In commencing the study, an attempt was made to optimize the launch date, 
along with several other parameters, for fixed flight time and final ecliptic inclina- 
tion. Considerable difficulty in convergence was encountered, however, until the 
launch date was fixed. A sequence of solutions over a range of fixed launch dates 
then indicated that the final mass was very insensitive to the launch date, varying 
only a few hundredths of a kilogram over a launch date range of 2-3 weeks. The 
convergence problem was caused by the fact that the time transversality was not a 
monotonic function of launch date. This situation causes the iterator to "hang-up 1 ' 
at the local extremum in the transversality condition, leading to a singularity in 
the boundary value problem. By forcing the solution away from the singular point 
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to the locality of the true optimum, convergence to the optimum launch date can be 
achieved. But since the improvement in final mass was negligible the additional 
effort required to obtain the overall optimum was felt to be unwarranted. Conse- 
quently, the launch date for all cases generated in this study was fixed at April 21, 
1979, which is about one month after the vernal equinox passage. 

Numerical difficulty was also encountered as a result of an attempt to drive 
the final state to correspond to a circular orbit of exactly 1 AU. A ground rule of 
the study was that at distances below 1 AU, the arrays are to be tilted such that 
the power generated equals that developed at 1 AU. The simulation of this effect 
results in a discontinuity in slope of the power factor y at the distance of 1 AU. 
The numerical problem arose as the solution neared convergence when neighboring 
trajectories would terminate on opposite sides of the 1 AU threshold. In effect, 
we were left with a problem not unlike that resulting from the infinitesimal thrust/ 
coast phases discussed in the preceding section in which the partial derivative 
matrices differed on opposite sides of the threshold. The problem was alleviated 
by choosing a final circular orbit radius of 1,001 AU which is sufficiently far 
from the point of discontinuity to eliminate the convergence difficulty, yet near 
enough to the desired radius that the performance results are valid. 

The other ground rules of the study were as follows: 

(1) a housekeeping power of 650 watts is to be provided; 

(2) total array power is 21 kw, leaving a reference power of 20.35 kw; 

(3) specific mass of the array is 15 kg/kw; 

(4) specific mass of the rest of the propulsion system is 15 kg/kw; 

(5) propellant tankage factor is 0.035; 

(6) total propulsion system efficiency is 0.63; 

(7) specific impulse is 3000 seconds. 

(8) optimize thrust direction, switch points, excess speed, geocentric de- 
clination, burn time and travel angle. 
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The results of the trajectory study are presented in the following table. 

The data are tabulated for the four flight times t^ as a function of the final ecliptic 

inclination i over the range of 45 to 60 degrees in increments of three degrees. 

The first several columns, including the initial adjoint variables X , X , X , 

x y z 
o J o o 

X* , X° and X' , the reference thrust acceleration g, the hyperbolic excess 
x y z 

o o o 

speed v , and the geocentric asymptote declination 6, comprise the independent 

parameters of the boundary value problem. These parameters are necessary to 

reproduce the case with HILTOP. The remaining columns contain, in order, the 

initial thrust offset angle a, the change in ecliptic longitude over the mission AX, 

burn time t^, initial spacecraft mass m^, final spacecraft mass m^, maximum 

solar distance encountered r , and minimum solar distance encountered r . . 

max min 
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EXTRA-ECLIPTIC MISSION PERFORMANCE STUDY 
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-.5391 

4.1913 

6936,9 

-39.684 
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-.7258 
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-6. 7294 
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